﻿using System;
using System.Collections;
using System.Collections.Generic;

namespace Mclib
{
	/// <summary>
	/// A class for implementing a 1D histogram.
	/// </summary>
	public class Histogram
	{
		/// <summary>
		/// Defines the bin boundaries.
		/// if BinSizesAreUniform, this variable is null.
		/// Lowest element is always -Inf.  Biggest element is always +Inf.
		/// Boundaries.Length = 1+Counts.Length.
		/// </summary>
		public readonly double[] Boundaries;
		/// <summary>
		/// Defines the bin centers.  The center of the first and last bins are defined as D/2 away from the adjacent bin boundary
		/// where D is the distance between the boundaries of the most adjacent pair of bins.
		/// </summary>
		public readonly double[] Centers;
		/// <summary>
		/// True if each bin is the same size.  False otherwise.
		/// </summary>
		public readonly bool BinSizesAreUniform;
		/// <summary>
		/// The lower boundary of the low bin (except for the first bin which is -Inf).
		/// </summary>
		public readonly double LowBinMin;
		/// <summary>
		/// The upper boundary of the high bin (except for the last bin which is +Inf).
		/// </summary>
		public readonly double HighBinMax;
		/// <summary>
		/// Defines the number of items in each bin.  An item "x" is in bin "iBin"
		/// if x is greater or equal to Boundaries[iBin] and less than Boundaries[iBin+1].
		/// </summary>
		public int[] Counts;
		private double binRange,binSpan,hBinSpan,invBinSpan;
		private int ibMax;
		/// <summary>
		/// Constructs an instances with nBins bins where BinSizesAreUniform will be true.
		/// The first and last bins catch anything outside the range [LowBinMin,HighBinMax).
		/// </summary>
		/// <param name="LowBinMin">The lower boundary of the low bin (except for the first bin which is -Inf).</param>
		/// <param name="HighBinMax">The upper boundary of the high bin (except for the last bin which is +Inf).</param>
		/// <param name="nBins">The number of bins including the first and last bin which catch evereything outside the range [LowBinMin,HighBinMax).</param>
		public Histogram(double LowBinMin,double HighBinMax,int nBins)
		{
			//	The static bin size
			this.BinSizesAreUniform = true;
			this.LowBinMin = LowBinMin;
			this.HighBinMax = HighBinMax;
			//	The range of x-values
			binRange = HighBinMax-LowBinMin;
			//	Each bin must span this range.
			binSpan = binRange/(double)(nBins-2);
			hBinSpan = 0.5*binSpan; // Half the range.
			invBinSpan = 1.0/binSpan;  // Multiplication is faster than division.
			//	Initialize counter
			this.Counts = new int[nBins];
			ibMax = nBins-1;
			//	Initialize the bin centers.
			Centers = new double[nBins];
			for (int i=0; i<nBins; i++)
				Centers[i] = LowBinMin+(double)i*binSpan-hBinSpan;
		}
		/// <summary>
		/// Registers the samples by incrementing the corresponding histogram bin counts.
		/// </summary>
		/// <param name="x">The data to be registered.</param>
		public void RegisterSamples(double[] x)
		{
			int i,iBin;
			if (BinSizesAreUniform)
			{
				for (i=0; i<x.Length; i++)
				{
					iBin = 1+(int)((x[i]-LowBinMin)*invBinSpan);
					if (iBin<=0)
						Counts[0]++;
					else if (iBin>=ibMax)
						Counts[ibMax]++;
					else
						Counts[iBin]++;
				}
			}
			else
			{
				throw new NotImplementedException("TO DO");
			}
		}
		/// <summary>
		/// Registers the samples taken from a specified column of x by incrementing the corresponding histogram bin counts.
		/// </summary>
		/// <param name="x">The data to be registered.</param>
		/// <param name="iCol">The 0-based dcolumn indexer</param>
		public void RegisterSamples(double[,] x, int iCol)
		{
			int i,iBin,len;
			len = x.GetLength(0);
			if (BinSizesAreUniform)
			{
				for (i=0; i<len; i++)
				{
					iBin = 1+(int)((x[i,iCol]-LowBinMin)*invBinSpan);
					if (iBin<=0)
						Counts[0]++;
					else if (iBin>=ibMax)
						Counts[ibMax]++;
					else
						Counts[iBin]++;
				}
			}
			else
			{
				throw new NotImplementedException("TO DO");
			}
		}
		/// <summary>
		/// Computes a total and a max for this.Counts.
		/// </summary>
		/// <param name="includeTails">Include the first and last bins that catch all samples outside [LowBinMin, HighBinMax).</param>
		/// <param name="total">The total number of samples.</param>
		/// <param name="max">The number of samples in the largest bin.</param>
		public void Summarize(bool includeTails,out int total,out int max)
		{
			total = 0;
			max = 0;
			int i,c;
			if (includeTails)
				for (i=0; i<Counts.Length; i++)
				{
					c = Counts[i];
					total += c;
					if (c > max)
						max = c;
				}
			else
				for (i=1; i<ibMax; i++)
				{
					c = Counts[i];
					total += c;
					if (c > max)
						max = c;
				}
		}
		/// <summary>
		/// Creates a density plot.  The maximum entry is 1.0.
		/// </summary>
		/// <param name="includeTails">Include the first and last bins that catch all samples outside [LowBinMin, HighBinMax).</param>
		/// <param name="max">The number of items in the largest bin (output from this.Summarize)</param>
		/// <returns>A density plot.  if (includeTails==true), the length is equal to Counts.Length, otherwise Counts.Length-2.</returns>
		public double[] Density(bool includeTails,int max)
		{
			double invMax = 1.0/(double)max; // multiplication is faster than division.
			int i;
			if (includeTails)
			{
				double[] output = new double[Counts.Length];
				for (i=0; i<output.Length; i++)
					output[i] = invMax*(double)Counts[i];
				return output;
			}
			else
			{
				double[] output = new double[Counts.Length-2];
				for (i=0; i<output.Length; i++)
					output[i] = invMax*(double)Counts[i+1];
				return output;
			}
		}
	}
	/// <summary>
	/// The class for implementing a 2D histogram.
	/// </summary>
	public class Histogram2
	{
		/// <summary>
		/// Defines the bin boundaries.  If BinSizesAreUniform, this variable is null.
		/// Boundaries[0] gives the bin boundaries for the y-axis.
		/// Boundaries[1] gives the bin boundaries for the x-axis.
		/// Lowest element is always -Inf.  Highest element is always +Inf.
		/// Boundaries.Length = 1+Counts.Length.
		/// </summary>
		public readonly double[][] Boundaries;
		/// <summary>
		/// Defines the bin centers.  The center of the first and last bins are defined as D/2 away from the adjacent bin boundary
		/// where D is the distance between the boundaries of the most adjacent pair of bins.
		/// Centers[0] gives the centers of every bin for the y-axis.
		/// Centers[1] gives the centers of every bin for the x-axis.
		/// </summary>
		public readonly double[][] Centers;
		/// <summary>
		/// True if each bin is the same size.  False otherwise.
		/// </summary>
		public readonly bool BinSizesAreUniform;
		/// <summary>
		/// The lower boundary of the low bin (except for the first bin which is -Inf).
		/// </summary>
		public readonly double[] LowBinMin;
		/// <summary>
		/// The upper boundary of the high bin (except for the last bin which is +Inf).
		/// </summary>
		public readonly double[] HighBinMax;
		/// <summary>
		/// Defines the number of items in each bin.  An item "x" is in bin "iBin"
		/// if x is greater or equal to Boundaries[iBin] and less than Boundaries[iBin+1].
		/// </summary>
		public int[,] Counts;
		private double[] binRange=new double[2],binSpan=new double[2],hBinSpan=new double[2],invBinSpan=new double[2];
		private int[] ibMax=new int[2];
		/// <summary>
		/// Constructs an instance where this.BinSizesAreUniform will be true.
		/// </summary>
		/// <param name="xLowBinMin">The lower boundary of the low bin for x (except for the first bin which is -Inf).</param>
		/// <param name="xHighBinMax">The upper boundary of the high bin for x (except for the last bin which is +Inf).</param>
		/// <param name="yLowBinMin">The lower boundary of the low bin for y (except for the first bin which is -Inf).</param>
		/// <param name="yHighBinMax">The upper boundary of the high bin for y (except for the last bin which is +Inf).</param>
		/// <param name="nxBins">The number of bins by the x-axis including the first and last bins which catch everything 
		/// out of the range specified by xLowBinMin and xHighBinMax.  The total number of bins will be nxBins*nyBins.</param>
		/// <param name="nyBins">The number of bins by the y-axisincluding the first and last bins which catch everything 
		/// out of the range specified by yLowBinMin and yHighBinMax.  The total number of bins will be nxBins*nyBins.</param>
		public Histogram2(double xLowBinMin,double xHighBinMax,double yLowBinMin,double yHighBinMax,int nxBins,int nyBins)
		{
			//	Set the basic parameters that define this 2D histogram.
			this.BinSizesAreUniform = true;
			this.LowBinMin = new double[] { xLowBinMin,yLowBinMin };
			this.HighBinMax = new double[] { xHighBinMax,yHighBinMax };
			//	Initialize counter
			this.Counts = new int[nyBins,nxBins];
			//	The range of values
			binRange[0] = HighBinMax[0]-LowBinMin[0];
			binRange[1] = HighBinMax[1]-LowBinMin[1];
			//	Each bin must span this range.
			binSpan[0] = binRange[0]/(double)(nyBins-2.0);
			binSpan[1] = binRange[1]/(double)(nyBins-2.0);
			// Half the range.
			hBinSpan[0] = 0.5*binSpan[0];
			hBinSpan[1] = 0.5*binSpan[1];
			// Multiplication is faster than division.
			invBinSpan[0] = 1.0/binSpan[0];
			invBinSpan[1] = 1.0/binSpan[1];
			ibMax[0] = nyBins-1;
			ibMax[1] = nxBins-1;
			//	Initialize the bin centers.
			int i,d,len;
			Centers = new double[2][];
			for (d=0; d<2; d++)
			{
				len = Counts.GetLength(d);
				Centers[d] = new double[len];
				for (i=0; i<len; i++)
					Centers[d][i] = LowBinMin[d]+(double)i*binSpan[d]-hBinSpan[d];
			}
		}
		/// <summary>
		/// Registers the samples taken from two specified columns of x by incrementing the corresponding histogram bin counts.
		/// </summary>
		/// <param name="x">The data to be registered.</param>
		/// <param name="iCol">The 0-based dcolumn indexer corresponding to the y-axis (row) of the histogram bin.</param>
		/// <param name="jCol">The 0-based dcolumn indexer corresponding to the x-axis (column) of the histogram bin.</param>
		public void RegisterSamples(double[,] x,int iCol,int jCol)
		{
			int i,iBin,jBin,len;
			len = x.GetLength(0);
			if (BinSizesAreUniform)
			{
				for (i=0; i<len; i++)
				{
					iBin = 1+(int)((x[i,iCol]-LowBinMin[0])*invBinSpan[0]);
					jBin = 1+(int)((x[i,jCol]-LowBinMin[1])*invBinSpan[1]);
					if (iBin<=0)
						iBin=0;
					else if (iBin>=ibMax[0])
						iBin = ibMax[0];
					if (jBin<=0)
						jBin=0;
					else if (jBin>=ibMax[1])
						jBin = ibMax[1];
					Counts[iBin,jBin]++;
				}
			}
			else
			{
				throw new NotImplementedException("TO DO");
			}
		}
		/// <summary>
		/// Computes a total and a max for this.Counts.
		/// </summary>
		/// <param name="includeTails">Include the first and last bins that catch all samples outside [LowBinMin, HighBinMax).</param>
		/// <param name="total">The total number of samples.</param>
		/// <param name="max">The number of samples in the largest bin.</param>
		public void Summarize(bool includeTails,out int total,out int max)
		{
			total = 0;
			max = 0;
			int iLen,jLen,iStart;
			if (includeTails)
			{
				iLen = Counts.GetLength(0);
				jLen = Counts.GetLength(1);
				iStart = 0;
			}
			else
			{
				iLen = ibMax[0];
				jLen = ibMax[1];
				iStart = 1;
			}
			int i,j,c;
			for (i=iStart; i<iLen; i++)
				for (j=iStart; j<jLen; j++)
				{
					c = this.Counts[i,j];
					if (c>max)
						max = c;
					total += c;
				}
		}
		/// <summary>
		/// Creates a density plot.  The maximum entry is 1.0.
		/// </summary>
		/// <param name="includeTails">Include the first and last bins that catch all samples outside [LowBinMin, HighBinMax).</param>
		/// <param name="max">The number of items in the largest bin (output from this.Summarize)</param>
		/// <returns>Density normalized so the maximum bin value is 1.0.</returns>
		public double[,] Density(bool includeTails,int max)
		{
			double invMax = 1.0/(double)max; // multiplication is faster than division.
			int i,j;
			int iLen,jLen,iStart;
			if (includeTails)
			{
				iLen = Counts.GetLength(0);
				jLen = Counts.GetLength(1);
				iStart = 0;
			}
			else
			{
				iLen = ibMax[0];
				jLen = ibMax[1];
				iStart = 1;
			}
			double[,] output = new double[iLen-iStart,jLen-iStart];
			for (i=iStart; i<iLen; i++)
				for (j=iStart; j<jLen; j++)
					output[i,j] = invMax*(double)Counts[i,j];
			return output;
		}
	}
}